#### Replication script for analyses in "Trade and Redistribution: Trade Politics and the Origins of Redistributive Taxation"

# log to file: 
# sink("psrm_replication_log.Rout")

### Preliminaries

## Set an appropriate working directory
# setwd(YOUR DIRECTORY) ## this should be where the data files "trade_redist_macro_tscs.RData" and "trade_redist_brit_constit.RData" are saved

## Load (install if necessary) required packages:
library(stargazer)
library(compactr)
library(clusterSEs)
library(multiwayvcov)
library(lmtest)
library(erer)
library(MASS)
library(foreign)
library(mvtnorm)

set.seed(12345)

## define inverse logit function

inv.logit <- function(p) {
  return(exp(p)/(1+exp(p)))
}
  
#### Macro Comparative Analyses: Tables 1 & 2; figures 3 & 4

####### Load data
load("trade_redist_macro_tscs.RData") # data is called macro


####### Table 1

market.empty = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq 
                  + as.factor(country_name) + as.factor(year), data = macro)

market.preferred = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                      + l1.e2 + l1.trade + l1.rgdppc + votetax2
                      + as.factor(country_name) + as.factor(year), data = macro)

market.full = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                 + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.taxgdp
                 + as.factor(country_name) + as.factor(year), data = macro)

market.no.int.preferred = lm(market ~ l1.tradeallyadv + l1rural.ineq
                             + l1.e2 + l1.trade + l1.rgdppc + votetax2
                             + as.factor(country_name) + as.factor(year), data = macro)

direct.empty = lm(direct~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq 
                  + as.factor(country_name) + as.factor(year), data = macro)

direct.preferred = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                      + l1.e2 + l1.trade + l1.rgdppc + votetax2
                      + as.factor(country_name) + as.factor(year), data = macro)

direct.no.int.preferred = lm(direct ~ l1.tradeallyadv + l1rural.ineq
                             + l1.e2 + l1.trade + l1.rgdppc + votetax2
                             + as.factor(country_name) + as.factor(year), data = macro)


direct.full = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                 + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.taxgdp
                 + as.factor(country_name) + as.factor(year), data = macro)

stargazer(market.empty, market.preferred, market.full, direct.empty, direct.preferred,direct.full, omit = c("year", "country", "Constant")
          , type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "taxgdp")
          , dep.var.labels = c("Domestic Indirect Share ", "Direct Share")
          , covariate.labels = c("Labour trade advantage (LTA)", "Inequality", "LTA : inequality", "GDP per capita", "Vote-tax link","Economic franchise",  "Trade",  "Tax/GDP")
          , style = "ajps", omit.stat = c("ser", "f"), add.lines = list(c("Country fixed effects", rep("Y", 6)), c("Year fixed effects", rep("Y", 6))), digits = 2
          , table.placement = "h!", title = "The determinants of revenue shares in total taxation, 1870-1913", label = "tab:main"
)

##### Table 2 -- top tax rates
ytax.empty = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                + as.factor(country_name) + as.factor(year), data = macro)

ytax.preferred = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                    + l1.e2 + l1.trade + l1.rgdppc + votetax2
                    + as.factor(country_name) + as.factor(year), data = macro)

ytax.preferred.noint = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq 
                          + l1.e2 + l1.trade + l1.rgdppc + votetax2
                          + as.factor(country_name) + as.factor(year), data = macro)

ytax.full = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
               + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.taxgdp
               + as.factor(country_name) + as.factor(year), data = macro)

htax.empty = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                + as.factor(country_name) + as.factor(year), data = macro)

htax.preferred = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                    + l1.e2 + l1.trade + l1.rgdppc + votetax2
                    + as.factor(country_name) + as.factor(year), data = macro)

htax.full = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
               + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.taxgdp
               + as.factor(country_name) + as.factor(year), data = macro)

stargazer(ytax.empty, ytax.preferred, ytax.full, htax.empty, htax.preferred, htax.full, omit = c("year", "country", "Constant")
          , type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "taxgdp")
          , dep.var.labels = c("Top Income Tax Rate ", "Top Inheritance Tax Rate")
          , covariate.labels = c("Labour trade advantage (LTA)", "Inequality", "LTA : inequality", "GDP per capita", "Vote-tax link","Economic franchise",  "Trade",  "Tax/GDP")
          , style = "ajps", omit.stat = c("ser", "f"), add.lines = list(c("Country fixed effects", rep("Y", 6)), c("Year fixed effects", rep("Y", 6))), digits = 2
          , table.placement = "p", title = "The determinants of top direct tax rates, 1870-1913", label = "tab:topRates"
)

### Figure 3
## figure 3a
beta.hat <- coef(market.preferred)
# pull out the covariance matrix
cov <- vcov(market.preferred)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(market.preferred$model$l1rural.ineq), max(market.preferred$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx <- beta.hat["l1.tradeallyadv"] + beta.hat["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx <- sqrt(cov["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

# plot the ME using compactr
#par(mfrow = c(1,1))
par(mar = c(4,5,2,2)+0.4)
plot(density(market.preferred$model$l1rural.ineq), type = "l",  xlim = mm(market.preferred$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(market.preferred$model$l1rural.ineq)

mnri = min(market.preferred$model$l1rural.ineq, na.rm = T)
mxri = max(market.preferred$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx, lwd = 3, type = "l", xlim = mm(market.preferred$model$l1rural.ineq), ylim = mm(c(upr, lwr)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"revenue share", partialdiff*"LTA")))
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr, lty = 3)
lines(z0, upr, lty = 3)

dev.print(device = pdf, file = "figure_3a.pdf", height = 6, width = 6)

## figure 3b
beta.hat <- coef(direct.preferred)
# pull out the covariance matrix
cov <- vcov(direct.preferred)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(direct.preferred$model$l1rural.ineq), max(direct.preferred$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx <- beta.hat["l1.tradeallyadv"] + beta.hat["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx <- sqrt(cov["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

# plot 
par(mar = c(4,5,2,2)+0.4)
plot(density(direct.preferred$model$l1rural.ineq), type = "l",  xlim = mm(direct.preferred$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(direct.preferred$model$l1rural.ineq)

mnri = min(direct.preferred$model$l1rural.ineq, na.rm = T)
mxri = max(direct.preferred$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx, lwd = 3, type = "l", xlim = mm(direct.preferred$model$l1rural.ineq), ylim = mm(c(upr, lwr)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"revenue share", partialdiff*"LTA")))
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr, lty = 3)
lines(z0, upr, lty = 3)


dev.print(device = pdf, file = "figure_3b.pdf", height = 6, width = 6)

###### Figure 4
###### Figure 4a:
beta.hat <- coef(ytax.preferred)
# pull out the covariance matrix
cov <- vcov(ytax.preferred)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(ytax.preferred$model$l1rural.ineq), max(ytax.preferred$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx <- beta.hat["l1.tradeallyadv"] + beta.hat["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx <- sqrt(cov["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

# plot the ME using compactr

#par(mfrow = c(1,1))
par(mar = c(4,5,2,2)+0.4)
plot(density(ytax.preferred$model$l1rural.ineq), type = "l",  xlim = mm(ytax.preferred$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(ytax.preferred$model$l1rural.ineq)

mnri = min(ytax.preferred$model$l1rural.ineq, na.rm = T)
mxri = max(ytax.preferred$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx, lwd = 3, type = "l", xlim = mm(ytax.preferred$model$l1rural.ineq), ylim = mm(c(upr, lwr)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"top rate", partialdiff*"LTA")))
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr, lty = 3)
lines(z0, upr, lty = 3)

dev.print(device = pdf, file = "figure_4a.pdf", height = 6, width = 6)

##### Figure 4b
beta.hat <- coef(htax.preferred)
# pull out the covariance matrix
cov <- vcov(htax.preferred)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(htax.preferred$model$l1rural.ineq), max(htax.preferred$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx <- beta.hat["l1.tradeallyadv"] + beta.hat["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx <- sqrt(cov["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

par(mar = c(4,5,2,2)+0.4)
plot(density(htax.preferred$model$l1rural.ineq), type = "l",  xlim = mm(htax.preferred$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(htax.preferred$model$l1rural.ineq)

mnri = min(htax.preferred$model$l1rural.ineq, na.rm = T)
mxri = max(htax.preferred$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx, lwd = 3, type = "l", xlim = mm(htax.preferred$model$l1rural.ineq), ylim = mm(c(upr, lwr)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"top rate", partialdiff*"LTA")))
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr, lty = 3)
lines(z0, upr, lty = 3)


dev.print(device = pdf, file = "figure_4b.pdf", height = 6, width = 6)

############## British Data
load("trade_redist_brit_constit.RData") ## data is called britcon

## Table 3 and Figure 5: ease of Labour entry

### models

ease1ft = polr(labelease ~ trade.gen.f + contest1900 + cvs1900.2+ cvs1900.22, method = "logistic", data = britcon)

ease2ft = polr(labelease ~ trade.gen.f + contest1900 + cvs1900.2 + cvs1900.22+ scotland, method = "logistic", data = britcon)

ease3ft = polr(labelease ~ trade.gen.f + lb.socec + contest1900 + cvs1900.2+ cvs1900.22 + scotland, method = "logistic", data = britcon, Hess = T)

ease4ft = polr(labelease ~ trade.gen.f + lb.socec + double + d.type + contest1900 + cvs1900.2 + cvs1900.22+ scotland, method = "logistic", data = britcon, Hess = T)

ease5ft = polr(labelease ~ trade.gen.f + lb.socec + double + d.type + industrial.f+ contest1900 + cvs1900.2+ cvs1900.22 + scotland, method = "logistic", data = britcon, Hess = T)

## print table to console (see Table 3)
stargazer(ease1ft, ease2ft, ease3ft, ease4ft, ease5ft
          , style = "ajps"
          , title="Ease of labor election in 1906. Model 4 provides estimates for figure~5"
          , label = "tab:labease1"
          , order = c("trade.gen.f1", "trade.gen.f0", "double", "typec", "socecB", "socecC", "industrial.f1", "industrial.f2","contest", "cvs", "scotland")
          , covariate.labels = c("Free trade", "Neutral trade", "Double member", "County", "Mixed class", "Working class", "Industrial", "Part-industrial","Election 1900: Liberal", "Election 1900: Conservative", "Election 1900: Neither", "Cons. vote 1900", "Cons. vote 1900 squared","Scotland")
          , notes.append = T
          , dep.var.labels = "Dependent variable: ordered factor."
          , notes = c("Variables beginning `Election' indicate constituencies contested by only one, or neither, of the major parties -- the omitted category contains Lib-Con contested constituencies.")
          , type = "text"
          , digits = 2
)

## Create marginal effects for Figure
ease1ftmfx = ocME(ease1ft)
ease2ftmfx = ocME(ease2ft)
ease3ftmfx = ocME(ease3ft)
ease4ftmfx = ocME(ease4ft)
ease5ftmfx = ocME(ease5ft)

locbarsft = barplot(ease4ftmfx$out$ME.0[,1], horiz = T, beside = T, las = 2)
lowerft0 = ease4ftmfx$out$ME.0[,1] - 1.96*ease4ftmfx$out$ME.0[,2]
upperft0 = ease4ftmfx$out$ME.0[,1] + 1.96*ease4ftmfx$out$ME.0[,2]
lowerft1 = ease4ftmfx$out$ME.1[,1] - 1.96*ease4ftmfx$out$ME.1[,2]
upperft1 = ease4ftmfx$out$ME.1[,1] + 1.96*ease4ftmfx$out$ME.1[,2]
lowerft2 = ease4ftmfx$out$ME.2[,1] - 1.96*ease4ftmfx$out$ME.2[,2]
upperft2 = ease4ftmfx$out$ME.2[,1] + 1.96*ease4ftmfx$out$ME.2[,2]

plotorder4ft = c(12, 11, 10, 9, 8, 7, 4, 3, 6, 5, 1, 2)
plotnames4ft = c("Scotland", "Cons. vote 1900^2", "Cons. vote 1900", "Election 1900: Neither", "Election 1900: Cons.", "Election 1900: Liberal", "Working class", "Mixed class", "County", "Double-member", "Neutral trade","Free trade")
par(mar = c(2.5, 0.7, 1, 0.8)+0.2, mfrow = c(1,3), oma = c(3, 10.5, 3, 0))

## Three-panel plot, Figure 5
plot(ease4ftmfx$out$ME.0[plotorder4ft,1], locbarsft, pch = 18, xlim = c(-0.6, 0.4), yaxt = "n", xlab = "", ylab = "", main = "")
abline(v = 0, col = "grey")
segments(x0 = lowerft0[plotorder4ft], x1 = upperft0[plotorder4ft], y0 = locbarsft)
#mtext(text = rownames(ease4mfx$out$ME.0)[plotorder4ft], at = locbars, line = 1, side =2, las = 2)
mtext(text = plotnames4ft, at = locbarsft, line = 1, side =2, las = 2, cex = 0.8)
mtext("No Labor \n Candidate", side = 3, line = 1, font = 3, cex = 0.8)

plot(ease4ftmfx$out$ME.1[plotorder4ft,1], locbarsft, pch = 18, xlim = c(-0.6, 0.4), yaxt = "n", xlab = "", ylab = "")
abline(v = 0, col = "grey")
segments(x0 = lowerft1[plotorder4ft], x1 = upperft1[plotorder4ft], y0 = locbarsft)
mtext("Labor Candidate, \n Opposed", side = 3, line = 1, font = 3, cex = 0.8)
mtext("Marginal effect of unit change \n Calculated at mean of independent variables", side = 1, outer = T, line = 1.5, cex = 0.8)

plot(ease4ftmfx$out$ME.2[plotorder4ft,1], locbarsft, pch = 18, xlim = c(-0.6, 0.4), yaxt = "n", xlab = "", ylab = "")
abline(v = 0, col = "grey")
segments(x0 = lowerft2[plotorder4ft], x1 = upperft2[plotorder4ft], y0 = locbarsft)
mtext("Labor Candidate, \n Unopposed", side = 3, line = 1, font = 3, cex = 0.8)

dev.print(file = "figure5.pdf",device = pdf, width = 7.5, height = 4)

### constituencies by tariff reform and party (fig 6)
dev.off()
rti_party = prop.table(table(britcon$trade.gen, britcon$party2), margin = 2)[ , c(2,9,10,11)]

rti_party_abs = table(britcon$trade.gen, britcon$party2)[ , c(2,9,10,11)]

shadesOfGrey <- colorRampPalette(c("grey0", "grey100"))
threeGreys <- shadesOfGrey(3)
fourGreys <- shadesOfGrey(4)

par(mar = c(4,6,1,1)+0.2)

pfiii=barplot(rti_party_abs[ , c(4,2,1,3)], beside = T, horiz = T, las = 1, xlim = c(0, 250), col = threeGreys, names.arg = c("Liberal \n Unionist", "Labour", "Conservative/\n Unionist", "Liberal"))
mtext("Number of constituencies", side = 1, line = 2.5)

legend(50, 8, legend = c("Free Trade", "Neutral / No position", "Tariff reform"), fill = rev(threeGreys), bty = "n")

dev.print(device = pdf, file = "figure6.pdf", width = 7.5, height = 4.5)

## main table models: 
##  cluster errors by constituency type 

## 0. trade interests only
m0f.type = glm(aye ~  trade.gen.f, family=binomial(link = "logit"), data = britcon)
m0f.t.vcov<-cluster.vcov(m0f.type, britcon$ctype)
m0f.t.coef = coeftest(m0f.type, m0f.t.vcov)

## 1. pelling controls (all categorical)
m1f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec, family=binomial(link = "logit"), data = britcon)
m1f.t.vcov<-cluster.vcov(m1f.type, britcon$ctype)
m1f.t.coef = coeftest(m1f.type, m1f.t.vcov)

## 2. plus party
m2f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec + party2, family=binomial(link = "logit"), data = britcon, x=T)
m2f.t.vcov<-cluster.vcov(m2f.type, britcon$ctype)
m2f.t.coef = coeftest(m2f.type, m2f.t.vcov)

## 3. plus region 

m3f.type = glm(aye ~  trade.gen.f + military + industrial.f + d.type + lb.socec + region2, family=binomial(link = "logit"), data = britcon, x=T)
m3f.t.vcov<-cluster.vcov(m3f.type, britcon$ctype)
m3f.t.coef = coeftest(m3f.type, m3f.t.vcov)

# party and region controls
m4f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec + party2 + region2, family=binomial(link = "logit"), data = britcon, x=T)
m4f.t.vcov<-cluster.vcov(m4f.type, britcon$ctype)
m4f.t.coef =coeftest(m4f.type, m4f.t.vcov)

## 5. party-region fes
m5f.type = glm(aye ~ trade.gen.f + lb.socec + industrial.f + d.type + military + rp, data = britcon, family=binomial(link = "logit"))
m5f.t.vcov<-cluster.vcov(m5f.type, britcon$ctype, force_posdef = T)
m5f.t.coef=coeftest(m5f.type, m5f.t.vcov)

## 6. limited to split delegations, party and region
m6f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec + party2 + region2, data = britcon[which(britcon$rp == "Southwest-Liberal" | britcon$rp == "South East-Liberal" | britcon$rp == "South East-Conservative/Unionist" | britcon$rp == "North-Conservative/Unionist" | britcon$rp == "East Midlands-Liberal" | britcon$rp == "Central-Liberal"),], family=binomial(link = "logit"))
m6f.t.vcov<-cluster.vcov(m6f.type, britcon$ctype[which(britcon$rp == "Southwest-Liberal" | britcon$rp == "South East-Liberal" | britcon$rp == "South East-Conservative/Unionist" | britcon$rp == "North-Conservative/Unionist" | britcon$rp == "East Midlands-Liberal" | britcon$rp == "Central-Liberal")])
m6f.t.coef=coeftest(m6f.type, m6f.t.vcov)

## 7. limited to split delegations, no party, no region 
m7f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec, data = britcon[which(britcon$rp == "Southwest-Liberal" | britcon$rp == "South East-Liberal" | britcon$rp == "South East-Conservative/Unionist" | britcon$rp == "North-Conservative/Unionist" | britcon$rp == "East Midlands-Liberal" | britcon$rp == "Central-Liberal"),], family=binomial(link = "logit"))
m7f.t.vcov<-cluster.vcov(m7f.type, britcon$ctype[which(britcon$rp == "Southwest-Liberal" | britcon$rp == "South East-Liberal" | britcon$rp == "South East-Conservative/Unionist" | britcon$rp == "North-Conservative/Unionist" | britcon$rp == "East Midlands-Liberal" | britcon$rp == "Central-Liberal")])
m7f.t.coef=coeftest(m7f.type, m7f.t.vcov)

## Table 4:

stargazer(m0f.type, m1f.type, m3f.type, m2f.type, m4f.type, m7f.type, m6f.type
          ,  se = list(m0f.t.coef[,2], m1f.t.coef[,2], m3f.t.coef[,2], m2f.t.coef[,2], m4f.t.coef[,2], m7f.t.coef[,2], m6f.t.coef[,2])
          , omit = c("Constant", "region2", "party2", "rp")
          , style = "ajps"
          , order = c("trade.gen.f1", "trade.gen.f0", "lb.socec", "industrial.f")
          , digits = 2
          , covariate.labels = c("Trade: free trade", "Trade: neutral", "Mixed socio-economic", "Working class", "Part-industrial", "Industrial", "Military", "County (more rural)") 
          , dep.var.labels = "`Aye' vote"
          , notes.append = T
          , notes = c("Standard errors clustered by constituency type.", "Omitted categories are: protectionist trade, upper class socioceonomic, non-industrial.", "`Split delegation' models include only MPs from parties and regions which included both `aye' and `no' votes.")
          , label = "tab:ukmain"
          , title = "Determinants of MPs' votes on the the Third Reading of the 1909 Finance Bill."
          , add.lines = list(c("Includes region", "No", "No", "Yes", "No", "Yes", "No", "Yes"), c("Includes party", "No", "No", "No", "Yes", "Yes", "No", "Yes"))
          , type = "text"
) 

#### Figure 7

m3f.sim <- rmvnorm(1000, m3f.type$coef, m3f.t.vcov)

freetrade.3cx <- c(1, 0,1,0,0,0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0) # liberal, not military, non-industrial, upper class borough; free trade interests
protection.3cx <- c(1, 0,0,0,0,0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0)  # liberal, not military, non-industrial, upper class borough; protectionist trade interests

pp.3freetrade <- numeric(1000)
for(j in 1:1000) {
  pp.3freetrade[j] <- t(as.matrix(freetrade.3cx)) %*% m3f.sim[j,]
}

pe.3ft <- mean(inv.logit(pp.3freetrade))
lo.3ft <- quantile(inv.logit(pp.3freetrade), .025)
hi.3ft <- quantile(inv.logit(pp.3freetrade), .975)

pp.3protection <- numeric(1000)
for(j in 1:1000) {
  pp.3protection[j] <- t(as.matrix(protection.3cx)) %*% m3f.sim[j,]
}

pe.3pr <- mean(inv.logit(pp.3protection))
lo.3pr <- quantile(inv.logit(pp.3protection), .025)
hi.3pr <- quantile(inv.logit(pp.3protection), .975)

# for error bars
lo90.3ft <- quantile(inv.logit(pp.3freetrade), .05)
hi90.3ft <- quantile(inv.logit(pp.3freetrade), .95)
lo90.3pr <- quantile(inv.logit(pp.3protection), .05)
hi90.3pr <- quantile(inv.logit(pp.3protection), .95)

par(mar = c(4,6,1,1)+0.2)
brep = barplot(height = c(pe.3ft, pe.3pr), horiz = T, xlim = c(-0.05, 1.05), xlab = "Expected probability of `aye'", col = fourGreys[3], names.arg = c("Free trade", "Protectionist"), las = 1 , ylim = c(-0.5, 4), space = 0.5)

segments(lo.3ft, brep[1], hi.3ft, brep[1], lty=3)
segments(lo.3pr, brep[2], hi.3pr, brep[2], lty=3)

points(x = c(lo.3ft,  hi.3ft), y = c(brep[1],brep[1]), pch = 18)
points(x = c(lo.3pr,  hi.3pr), y = c(brep[2],brep[2]), pch = 18)

dev.print(device = pdf, file = "figure7.pdf", width = 5.5, height = 1.669875)

#sink()
